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Abstract 

This  work  used  atomistic  molecular  dynamics  simulations  to  predict  envi¬ 
ronmental  impact  of  six  energetic  materials,  2,4-dinitroanisole  (DNAN), 
N-methyl-p-nitroaniline  (MNA),  3,5-dinitropyrazole  (DNP),  3-nitro-i,2,4- 
triazol-5-one  (NTO),  i-methyl-2,4,5-trinitroimidazole  (MTNI)  and  1,3,5- 
triamino-2, 4, 6-trinitrobenzene  (TATB).  Molecular  models  developed  for 
these  compounds  were  used  to  determine  octanol-water  partition  coeffi¬ 
cient  (log  Kow)  and  Henry’s  law  constant  (log  H).  Log  Kow  was  predicted 
for  DNAN  and  MNA  to  within  ±0.1  log  units  of  experiment,  while  log  H 
was  predicted  to  within  ±1.0  log  units.  For  the  remaining  four  compounds, 
no  experimental  data  exist  for  comparison.  Predicted  log  Kow  and  log  H 
values  suggest  that  these  compounds  have  the  potential  to  cause  ground- 
water  contamination.  Depending  on  the  values  of  the  partition  coeffi¬ 
cients,  appropriate  treatment  methodologies  can  be  chosen  for  each  con¬ 
taminant  of  interest.  In  addition  to  partition  coefficients,  a  variety  of 
thermophysical  properties  were  predicted,  including  vapor-liquid  coexist¬ 
ence  curves,  critical  points,  vapor  pressure,  heats  of  vaporization,  crystal 
lattice  parameters,  and  solid  density.  The  crystal  density  and  lattice  pa¬ 
rameters  predicted  for  all  energetic  materials  were  in  close  agreement  with 
experimental  data.  Overall,  these  results  suggest  that  empirical  force 
fields,  combined  with  molecular  dynamics  simulations,  provide  an  accu¬ 
rate  methodology  for  predicting  relevant  descriptors  of  environmental  fate 
for  energetic  materials. 
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1  Introduction 

1.1  Background 

Over  the  past  two  decades,  there  has  been  considerable  work  in  the  devel¬ 
opment  of  insensitive  munitions  (IM),  which  exhibit  low  shock  sensitivity 
and  high  thermal  stability  due  to  the  increased  safety  and  environmental 
concerns  associated  with  the  traditionally  used  explosives.  Contamination 
of  the  environment  (ground  water,  soil,  and  sediment)  by  explosives  due  to 
military  activities  such  as  weapon  production  and  handling,  weapons  test¬ 
ing  and  training,  waste  discharge,  and  demilitarization  has  become  a  mul¬ 
ti-billion  dollar  problem.  Therefore,  development  of  methodologies  that 
can  be  used  to  predict  the  environmental  fate  of  a  particular  compound 
before  its  deployment  or  even  pre-synthesis  is  of  great  importance,  and 
could  lead  to  significant  long-term  cost  savings. 

The  fate  of  any  compound  in  soil,  water,  or  the  atmosphere  can  be  deter¬ 
mined  by  studying  the  interaction  between  the  compound  and  the  target 
medium.  These  interactions  are  described  in  part  by  partitioning  the  com¬ 
pound  of  interest  between  two  different  media,  which  is  represented  by 
various  partition  or  distribution  coefficients.  Two  key  partition  coefficients 
used  to  assess  a  compound’s  impact  on  air,  water,  and  organic  media  are 
octanol-water  partition  coefficients  and  Henry’s  law  constants  (air-water 
partition  coefficient).  The  Henry’s  law  constant  is  the  equilibrium  distribu¬ 
tion  of  a  species  between  gas  and  liquid.  For  dilute  aqueous  solutions, 
Henry’s  law  constant  is  the  ratio  of  the  solute’s  partial  pressure  and  its 
aqueous  concentration.  Higher  Henry’s  law  constant  means  higher  volatil¬ 
ity  and  lower  aqueous  solubility.  Octanol-water  partition  coefficient,  which 
is  the  ratio  of  the  concentration  of  a  neutral  chemical  species  in  octanol 
and  in  water  at  equilibrium,  is  a  measure  of  hydrophobicity/lipophilicity 
or  hydrophilicity  of  a  compound.  A  chemical  species  can  be  classified  as 
hydrophobic  (log  KoW  >6)  or  hydrophilic  (log  Kow  <o)  depending  on  the 
value  of  the  octanol-water  partition  coefficient. 

A  wide  variety  of  experimental  techniques  (Sangster  1997)  exist  to  meas¬ 
ure  partition  coefficients,  but  theoretical  methods  offer  a  few  benefits  over 
experiments  for  energetic  materials.  With  appropriate  computational 
methodologies,  it  is  possible  to  predict  the  behavior  of  a  material  in  the 
environment  before  it  has  been  synthesized,  allowing  for  a  prescreening  of 
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potential  candidate  molecules.  Such  a  prescreening  is  expected  to  lead  to 
cost  savings  by  reducing  the  pool  of  candidate  molecules  for  synthesis  and 
eliminating  the  need  for  environmental  remediation  near  manufacturing 
sites  and  test  ranges.  Moreover  the  hazardous  nature  and  long  experi¬ 
mental  time  scales  associated  with  the  development  and  testing  of  each 
compound  makes  computational  methods  an  alternative  for  experiments. 

Figure  l  shows  the  molecular  structures  of  the  six  energetic  materials.  In 
this  context,  molecular  models  or  “force  fields”  have  been  developed  for  all 
six  compounds  to  predict  octanol-water  partition  coefficient,  Henry’s  law 
constant,  vapor-liquid  coexistence  curves,  critical  parameters,  vapor  pres¬ 
sure,  boiling  point,  acentric  factor,  heats  of  vaporization,  lattice  parame¬ 
ters,  and  crystal  density.  The  motivation  for  this  study  comes  from  the 
work  of  James  W.  Gillett  who  proposed  a  comprehensive  prebiological 
screen  (Gillett  1983)  that  correlates  the  octanol-water  partition  coefficient 
and  the  Henry’s  law  constant  to  predict  which  materials  have  the  potential 
to  be  problematic  if  released  to  the  environment. 

Numerous  theoretical 
methods  have  been  re¬ 
ported  for  the  prediction 
of  partition  coefficients. 

Most  of  the  computational 
methods  are  based  on 
fragment/group  or  bond 
contribution  methodology 
(Leo  1993;  Klopman  et  al. 

1981;  Rekker  et  al.  1982; 

Suzuki  et  al.  1990;  Rekker 
et  al.  1979;  Broto  et  al. 

1984;  Ghose  et  al.  1986), 
which  use  molecular  de¬ 
scriptors  (topological,  topographical  and  quantum  chemical)  derived  from 
a  training  set  of  compounds.  Other  methods  include  Quantitative  Struc¬ 
ture-Activity  Relationship  (QSAR)  and  Quantitative  Structure-Property 
Relationship  (QSPR)  models  (Hansch  et  al.  1967;  Hansch  et  al.  1995; 
Bodor  et  al.  1989;  Kantola  et  al.  1991;  Famini  et  al.  1992;  Moriguchi  et  al. 
1992;  Ghasemi  et  al.  2007),  which  relate  the  molecular  structures  to  bio¬ 
logical  activity  or  other  physical  properties.  These  activities  or  properties 
are  expressed  as  a  function  of  the  partition  coefficients. 


1.  Molecular  structures  of  the  energetic  materials 
studied  in  this  work. 


NTO  MTNI  TATES 
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Although  these  methods  are  faster  than  experiments,  the  need  for  numer¬ 
ous  empirical  parameters  limits  their  predictive  capability  to  molecules 
with  strong  similarities  to  those  used  in  the  training  set.  Other  theoretical 
methods  offer  promise  as  alternatives  to  traditional  QSPR,  e.g.,  continuum 
solvent  methods  such  as  COSMO  (conductor  screening  module),  SM 
(Solvation  Model),  GB/SA  (Generalized  Born/Surface  Area  models)  (Still 
et  al.  1990;  Jean-Charles  et  al.  1991;  Hawkins  et  al.  1998;  Klamt  et  al. 

1998;  Giesen  et  al.  1996;  Best  et  al.  1997)  and  explicit  solvent  methods 
with  molecular  mechanics  force  fields  such  as  Monte  Carlo  (MC)  or  Mo¬ 
lecular  Dynamics  (MD)  coupled  with  free  energy  perturbation  (FEP) 
(Kollman  1993;  Straatsma  et  al.  1992).  For  flexible  systems,  the  preferred 
method  is  MC  or  MD  simulations,  since  these  methods  allow  averaging 
over  numerous  molecular  conformations.  This  work  discusses  the  applica¬ 
tion  of  MD  simulation  to  predict  the  pharmacokinetic  properties  of  ener¬ 
getic  materials. 

1.2  Objectives 

The  objective  of  this  research  was  to  use  molecular  simulations  to  deter¬ 
mine  a  variety  of  physical  properties  that  can  be  used  to  predict  the  envi¬ 
ronmental  fate  of  six  IM  compounds:  2,4-dinitroanisole  (DNAN),  N- 
methyl-p-nitroaniline  (MNA),  3,5-dinitropyrazole  (DNP),  3-nitro-i,2,4- 
triazol-5-one  (NTO),  i-methyl-2,4,5-trinitroimidazole  (MTNI)  and  1,3,5- 
triamino-2, 4, 6-trinitrobenzene  (TATB) . 

1.3  Approach 

This  work  accomplished  screening  by  locating  coordinates  for  compounds 
in  the  two  dimensional  mobility  and  multimedia  exposure  plot  (log  Kow  vs 
log  H),  which  is  divided  into  specific  regions,  each  characterized  by  a 
unique  ecotoxicologic  risk  or  concern  like  bioaccumulation,  ground  water 
pollution,  and  some  indirect  atmospheric  effects  like  ozone  depletion. 
Therefore  knowledge  of  both  the  quantities  identifies  potential  environ¬ 
mental  concerns  associated  with  each  material. 

1.4  Mode  of  technology  transfer 

This  report  will  be  made  accessible  through  the  World  Wide  Web  (WWW) 
at  URL:  http://libweb.erdc.usace.armv.mil 
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2  Conformational  Analysis 

It  is  believed  that  the  insensitivity  and  thermal  stability  of  these  IM  com¬ 
pounds  are  an  outcome  of  the  intramolecular  and  intermolecular  interac¬ 
tions.  These  explosives  derive  most  of  their  characteristics  from  the  nitro 
and  the  amino  functional  groups.  Therefore  the  prediction  of  rotational  bar¬ 
riers  offers  valuable  insight  into  the  strength  of  these  intramolecular  inter¬ 
actions.  The  internal  rotation  mechanism  around  the  C-N  bond  also  yields 
important  details  of  the  bond  rupture,  which  is  a  crucial  phenomenon  in  the 
decomposition  reaction  of  explosives.  These  data  are  also  required  for  the 
development  of  atomistic  force  fields  for  use  in  MD  simulations. 

2.1  DNAN  and  IVINA 

The  conformational  behavior  of  DNAN  and  MNA  was  analyzed  with 
Hartree-Fock  (HF),  Moller  Plesset  (MP2)  and  density  functional  theory 
(Levine  1991)  using  the  hybrid  B3LYP  functional  to  check  for  other  stable 
conformers  and  to  obtain  torsional  barriers.  The  6-3ig+(d,p)  basis  set  was 
used  for  all  calculations. 

2.1.1  Equilibrium  structures 

Figure  2  shows  a  schematic  of  DNAN  and  MNA  structures.  For  DNAN, 
MP2  and  B3LYP  calculations  predicted  a  co-planar  structure  with  the 
methoxy  group  and  the  p-nitro  group  in  plane  with  the  aromatic  ring, 
while  the  ortho-nitro  group  was  tilted  out  of  plane.  The  methoxy  group 
adopted  a  conformation  anti  to  the  o-nitro  group  to  avoid  steric  hindrance. 
In  contrast,  at  the  HF  level  of  theory,  the  methoxy  group  was  nearly  or¬ 
thogonal  to  the  plane  of  the  aromatic  ring,  with  the  p-nitro  group  in  plane 
with  the  ring,  and  the  o-nitro  group  out  of  plane  with  the  ring. 

The  minimum  energy  conformer  at  the  MP2  and  B3LYP  theories  is  the  one 
with  the  methoxy  group  co-planar  with  the  ring.  At  the  HF  theory,  a  per¬ 
pendicular  methoxy  group  conformation  was  found  to  be  the  minimum 
energy  conformation.  The  O1-C1-C2  angle  is  121.4  degrees  at  HF  and  de¬ 
creases  as  the  theory  level  increases  (HF>B3LYP>  MP2).  The  reverse 
(HF<B3LYP<MP2)  holds  true  for  the  O1-C1-C6  angle. 
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Figure  4.  Resonance  structures  for  MNA. 


For  MNA,  all  theories  predicted  a  planar  structure.  The  equilibrium  pa¬ 
rameters  agree  well  with  each  other  and  with  the  experimental  crystal 
structure  (Panunto  et  al.  1987).  The  substitution  of  a  methyl  group  in  place 
of  hydrogen  in  the  amine  group  and  conjugate  effects  among  the  strong 
electron-donor  amine  group  and  the  phenyl  ring  causes  the  nitrogen  to 
adopt  a  planar,  rather  than  pyramidal,  conformation.  The  equilibrium 
C-NH  bond  length  was  shorter  than  the  typical  equilibrium  C-N  single 
bond  of  1.45  A.  The  reason  for  shortening  of  the  C-NH  bond  is  due  to  the 
presence  of  some  double  bond  characteristics,  which  are  caused  by  the 
conjugate  effects  between  the  ring  and  the  amino  and  nitro  groups  (Figure 
4).  Similar  to  DNAN,  very  weak  intramolecular  hydrogen  bonding  (2.4  A) 
occurs  in  MNA  for  all  the  optimized  structures. 

2.1.2  Torsional  barriers 

Table  1  lists  the  predicted  barriers  to  rotation  in  DNAN  and  MNA.  Figures 
5,  6,  and  7,  respectively,  show  torsional  barriers  for  the  methoxy,  para- 
and  ortho-  nitro  groups  in  DNAN.  Figures  8  and  9  show  torsional  barriers 
for  amino  group  and  p-nitro  group  in  MNA.  For  all  the  functional  groups, 
B3LYP  and  HF  predicted  torsional  barriers  higher  than  the  MP2  barriers 
except  for  the  methylamine  group  where  MP2  predicted  a  lower  barrier 
than  B3LYP  and  HF.  This  is  due  to  the  electron  correlation  effects,  which 
are  significant  in  these  molecules. 


Table  1.  Rotational  barriers  in  kcal/mol  for  DNAN  and  MNA. 


Theory  Levels 

DNAN 

MNA 

Methoxy 

O-nitro 

P-nitro 

Amine 

P-nitro 

HF 

4.8 

1.2 

7.2 

6.9 

8.8 

B3LYP 

4.8 

1.0 

7.0 

9.9 

8.6 

MP2 

3.1 

1.7 

4.2 

6.1 

4.6 
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Figure  5.  Torsional  barrier  methoxy  group  (C-O-C-C)  in  DNAN; 
B3LYP  (black),  MP2  (red),  and  HF  (green). 


Figure  6.  Torsional  barrier  for  p-nitro  group  (ONCC)  in  DNAN; 
B3LYP  (black),  MP2  (red),  and  HF  (green). 


Figure  7.  Torsional  barrier  for  o-nitro  group  (ONCC)  in 
DNAN;  B3LYP  (black),  MP2  (red),  and  HF  (green). 
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Figure  8.  Torsional  barrier  for  methylamine  group  (C-N-C-C)  in 
MNA;  B3LYP  (black),  MP2  (red),  and  HF  (green). 


Figure  9.  Torsional  barrier  for  nitro  group  (O-N-C-C)  in 
MNA;  B3LYP  (black),  MP2(red),  and  HF  (green). 


The  scattered  torsional  curve  obtained  for  methoxy  torsion  is  a  result  of 
steric  crowding  effects  between  the  methoxy  and  the  bulky  ortho  substitu¬ 
ent  (nitro)  group  and  conjugation  between  lone  electron  pairs  on  the  oxy¬ 
gen  atom  and  the  aromatic  n  system.  The  ortho-nitro  group  resists  the  ro¬ 
tation  of  the  methoxy  group  and  eventually  tilts  when  the  methyl  group 
approaches  it  to  avoid  overlap.  This  is  evident  from  the  C3-C2-N1  angle, 
which  decreases  during  internal  rotation  of  methoxy  group,  which  indi¬ 
cates  the  drift  of  the  ortho-nitro  group,  and  which  is  manifested  in  the  plot 
of  rotational  barriers  with  multiple  local  minima  and  maxima. 

The  barrier  to  torsion  for  amino  group  in  MNA  shows  a  discontinuity  in 
energy  around  a  dihedral  angle  of  -50  and  150  degrees.  This  behavior  was 
also  observed  and  explained  for  methylamine  rotation  by  Birkett  et  al. 


ERDC/CERL  TR-10-26 


9 


(2002)  in  their  work  with  substituted  triazine  rings.  It  is  due  to  the  ability 
of  nitrogen  in  a  substituted  amine  group  to  be  both  planar  and  pyramidal; 
when  the  pyramidal  nature  changes,  there  is  a  significant  drop  in  energy. 

As  Figure  8  shows,  at  a  dihedral  angle  of  o  degrees,  MP2  predicts  non-zero 
energy  while  B3LYP  and  HF  theories  predict  o  degrees  as  the  lowest  ener¬ 
gy  conformer.  On  the  contrary,  the  equilibrium  structure  at  the  MP2  level 
has  minimum  energy  for  a  co-planar  methyl-amino  group  (dihedral  angle 
of  o  degrees). 

To  further  investigate  this  behavior,  MP2  calculations  were  run  with  a 
double  diffuse  function  (++)  and  larger  basis  set  (6-3iig+(d,p)),  but  both 
gave  similar  relative  energies.  Calculations  performed  with  Quadratic  Con¬ 
figuration  Interaction  Singles  Doubles  (QCISD)  theory  and  3-2ig  basis  set 
predict  the  o  degrees  dihedral  as  the  lowest  energy  conformer,  in  agree¬ 
ment  with  HF  and  B3LYP  results.  These  results  suggest  that,  for  molecules 
that  have  resonance  structures  such  as  MNA,  MP2  theory  may  give  erro¬ 
neous  results  for  the  lowest  energy  conformer. 

2.2  DNP  and  NTO 

2.2.1  Equilibrium  structures 

The  optimized  structures  of  DNP  and  NTO  at  HF,  B3LYP  and  MP2  levels 
of  theories  are  all  planar  with  respect  to  the  nitro  groups.  The  molecular 
parameters  for  optimized  NTO  at  all  three  levels  of  theories  agree  well 
with  each  other  and  the  crystal  structure  from  experiment  (Bolotina  et  al. 
2005).  No  experimental  structure  exists  for  DNP  to  make  a  comparison. 
Figure  10  shows  a  schematic  of  DNP  and  NTO. 


2.2.2  Torsional  barriers 

Figures  11, 12,  and  13  show  the  torsional  barriers  determined  for  the  nitro 
groups  in  both  the  compounds. 


Figure  10.  Schematic  of  DNP  (left)  and  NTO  (right). 
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Figure  11.  Torsional  barriers  for  nitro  group  (01-N1-C2-N2)  in 
DNP;  B3LYP  (black),  MP2  (red),  and  HF  (green). 


B3LYP  (black),  MP2  (red),  and  HF  (green). 


Figure  13.  Torsional  barriers  for  nitro  group  (01-N2-C1-N3)  in  NTO; 
B3LYP  (black),  MP2  (red),  and  HF  (green). 
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Table  2  lists  the  predicted  barriers  to  rotation  around  the  C-N  bond  in 
DNP  and  NTO.  For  DNP,  the  predicted  barrier  to  rotation  around  the  N4- 
C3  bond  is  higher  than  the  rotation  around  N1-C2  bond.  This  is  due  to  the 
location  of  the  nitro  group  in  each  case.  The  nitro  group  at  C3  is  located 
adjacent  to  the  amide  hydrogen  and  involves  hydrogen  bonding  whereas 
the  nitro  group  at  C2  has  no  amide  hydrogens  to  bond. 


Table  2.  Barriers  to  rotation  in  DNP  and  NTO. 


Theory  Levels 

DNP 

NTO 

01-N1-C2-N2 

03-N4-C3-N3 

01-N2-C1-N3 

HF 

5.5 

9.4 

6.6 

B3LYP 

4.8 

8.8 

7.6 

MP2 

3.2 

6.9 

6.0 
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3  Force  Field  Development 


Force  fields  can  be  conveniently  split  into  two  types  of  interactions,  bond¬ 
ed  and  non-bonded.  Bonded  interactions  account  for  the  conformational 
structure  of  the  molecule,  and  include  bond  stretching,  bond  bending,  and 
torsional  rotation  around  the  various  bonds.  Non-bonded  interactions  de¬ 
scribe  the  energetics  of  atom-atom  interactions  and  are  usually  described 
by  an  atom-atom  pair  interaction  potential. 

3.1  Non-bonded  interactions 


Non-bonded  interactions  between  atoms  in  each  molecule  were  represent¬ 
ed  with  a  standard  12-6  Lennard-Jones  potential  with  a  coulombic  term 
for  partial  charges: 


U(rf)  =  4s„. 
where: 


(  \ 

12 

f  \ 

6' 

a- 

a- 

y 

_ 

y 

r 

r 

l  y  J 

{  y  J 

Mi 


4  ™orv 


(1) 


nj  -  atom-atom  separation 
sij  =  LJ  well  depth 
an  —  LJ  diameter 
qi  =  partial  charge  on  atom  i 
qj  =  partial  charge  on  atom  j 
so  —  permittivity  of  vacuum. 


Cross  interaction  parameters  for  unlike  atoms  were  determined  through 
Lorentz-Berthelot  combining  rules  (Lorentz  1881;  Berthelot  1898): 


a..  = 

y 


s..  = 
y 


—  (0..  +  0 ..) 
2  !1  ji 


S..S  .. 
»  ]] 


(2) 

(3) 


Initial  estimates  of  the  partial  charges  for  each  molecule  were  determined 
through  a  CHELPG  (CHarges  from  Electrostatic  Potentials  using  a  Grid 
based  method)  analysis  by  fitting  to  a  electrostatic  potential  determined 
from  ab  initio  calculations  performed  at  the  HF/6-3ig+(d,p)  level  of  theo¬ 
ry  and  basis  set  with  Gaussian  03  (Gaussian  2003).  These  partial  charges 
were  rescaled  by  a  factor  of  0.94  to  improve  the  reproduction  of  experi¬ 
mentally  determined  octanol-water  partition  coefficients.  Two  different 
force  fields  were  developed  for  DNAN  and  MNA:  united-atom  (UA)  and 
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explicit  hydrogen  (EH).  In  the  UA  force  field,  all  hydrogens  bonded  to  car¬ 
bon  atoms  are  combined  with  carbon  to  form  a  single  interaction  site  (a 
pseudo-atom)  centered  on  the  nucleus  of  the  carbon  atom.  In  the  EH  force 
field,  all  atoms  are  modeled  explicitly,  with  their  interaction  sites  centered 
on  the  respective  atomic  nuclei.  The  EH  force  field  for  DNAN  and  MNA 
was  motivated  by  the  poor  performance  of  the  UA  force  field  in  the  predic¬ 
tion  of  crystal  lattice  parameters  and  solid  densities.  For  all  other  com¬ 
pounds,  only  an  EH  force  field  was  constructed.  In  this  report,  the  united- 
atom  force  fields  are  referred  to  as  DNAN-UA  and  MNA-UA  and  explicit 
hydrogen  force  fields  as  DNAN-EH  and  MNA-EH. 

For  the  united- atom  force  field,  Lennard- Jones  parameters  a  and  s  for 
each  interaction  site  were  transferred  from  analogous  compounds  previ¬ 
ously  parameterized  in  the  development  of  the  TraPPE-UA  force  field 
(Martin  et  al.  1998;  Wick  et  al.  2000;  Stubbs  et  al.  2004;  Wick  et  al. 

2005).  In  the  EH  version,  the  Lennard-Jones  parameters  for  the  nitro 
group  were  transferred  from  the  explicit  model  of  nitrobenzene  reported 
in  the  recent  work  by  Siepmann  and  co-workers  (Rai  et  al.  2008)  and  the 
rest  from  the  TraPPE  force  field  (Rai  et  al.  2007)  for  five-membered  rings. 
The  aromatic  ring  was  modeled  as  EH  wherever  necessary.  The  parame¬ 
ters  for  the  ring  were  transferred  from  explicit  model  of  benzene  (Rai  et  al. 
2008).  Tables  A1-A6  (pp  40-41 )  in  the  Appendix  A  to  this  report  list  the 
Lennard-Jones  parameters  and  partial  charges. 

3.2  Bonded  interactions 

A  harmonic  term  was  used  to  represent  interactions  due  to  bond  stretching: 

Ubond=K(r~ro)2  {4; 

where: 

kb  =  force  constant 
r  =  measured  bond  length 
r0  =  equilibrium  bond  length. 

Bond  angle  bending  is  also  represented  by  a  harmonic  potential: 

Ubad=K(6-6 0)2  (5; 

where: 

kd  =  force  constant 
9  =  measured  bond  length 
Oo  =  equilibrium  bond  angle. 


ERDC/CERL  TR-10-26 


14 


Barriers  to  rotation  about  various  dihedral  angles  were  controlled  through 
a  cosine  series  fit  to  ab  initio  calculations: 

Utors  =  XkJl  +  cos(n4>-f)]  (6 

where: 

k#  =  force  constant 
<f)  =  dihedral  angle 
n  =  multiplicity 
/  =  phase  angle. 
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4  Methodology  and  Simulation  Details 

4.1  Partition  coefficients 

4.1.1  Octanol-water  partition  coefficient 

The  octanol-water  partition  coefficient  (log  kow)  is  related  to  the  free  ener¬ 
gy  of  transfer  for  the  solute  between  water  and  water-saturated  octanol 
phase  by: 

AG  =  — 2.303 RT  log  kow 


where: 

R  =  universal  gas  constant 
T  =  temperature. 

Direct  calculation  of  free  energies  of 
transfer  between  water  and  octanol 
phases  is  possible  for  small  solutes 
(Martin  et  al.  1997),  but  is  extreme¬ 
ly  difficult  for  the  larger,  multi¬ 
functional  molecules  of  interest  in 
this  work.  Fortunately,  because  the 
Gibbs  free  energy  is  a  state  func¬ 
tion,  it  is  still  possible  to  calculate 
the  Gibbs  free  energy  of  transfer  of  phases  through  the  suitable  choice  of 
path.  In  this  work,  AG  was  computed  via  the  thermodynamic  path  where 
solute  A  is  slowly  transformed  to  solute  B  in  water  and  water-saturated 
octanol  (Figure  14). 

This  path  provides  a  means  for  calculating  the  relative  Gibbs  free  energy  of 
transfer,  which  is  defined  by: 

AAGTAB  =  AGTB-AGTA  =  AGTf  JA0  B)-AGiy,  XA&B) 

Tr  Tr  Tr  Tr(octjK  y  Tr(w)K  /o\ 


Figure  14.  Thermodynamic  cycle  used  to 
calculate  octanol-water  partition  coefficient. 
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where: 

AGn  =  free  energy  of  transformation. 
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The  relative  partition  coefficient  is  now  expressed  as: 


A  loqk 

^  ow 


-AA  GTrAB 

2.303 RT 


(9) 


Free  energy  differences  are  calculated  by  the  FEP  technique  (Kollman 
1993;  Straatsma  et  al.  1992)  combined  with  constant  pressure-tempera¬ 
ture  MD.  The  FEP  method  involves  slowly  transforming  solute  A  to  solute 
B  (either  A  or  B  is  the  compound  of  interest)  by  scaling  the  interaction  po¬ 
tential  through: 

U(X)  =  XUB  +  (1-X)UA 


where: 

A  =scaling  parameter  (value  between  o  and  1). 

The  FEP  method  allows  calculation  of  the  relative  Gibbs  free  energy  of 
transfer  AAG  from  which  the  relative  octanol-water  partition  coefficient 
(A  log  kow)  is  obtained.  The  absolute  partition  coefficient  of  target  molecule 
B  is  then  calculated  from  the  reference  molecule  A  from: 

logk  (B)  =  Alogk  +logk  (A) 

u  ow  u  ow  u  ow  (11) 


4.1.2  Henry’s  Law  constant 


The  Henry’s  law  constant  is  the  equilibrium  distribution  of  a  species  be¬ 
tween  gas  and  liquid.  For  dilute  aqueous  solutions,  it  is  the  ratio  of  the  so¬ 
lute’s  partial  pressure  and  its  aqueous  concentration.  The  Henry’s  law  con¬ 
stant,  expressed  in  terms  of  solvation  energy  of  a  solute  in  water,  is  given 
by  (Lin  et  al.  2002): 


AG 


*  sol 
i/W 


RT  In  10 


+  log. 


rtp: 

N , 


(12) 


where: 

AGi/w*S0  =  solvation  free  energy  of  species  i  in  solvent  water 
pw°  =  number  density  of  pure  water 
Na  =  Avagadro’s  number. 


The  solvation  free  energy  of  solute  i  in  water,  or  the  hydration  free  energy, 
is  the  free  energy  associated  with  the  transfer  of  solute  from  vacuum  to 
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water.  Similar  to  the  octanol-water  partition  coefficient,  a  thermodynamic 
path  is  constructed,  but  the  water-saturated  octanol  phase  is  replaced  by 
the  vacuum  phase.  Solute  i  is  transformed  to  j  in  both  water  and  vacuum. 
A  relative  Henry’s  law  constant  term  can  be  derived  using  Equation  to  and 
the  equation  for  the  Henry’s  law  constant  for  solute  j  given  by: 


log 


A  G*s°' 

_  j/w 

10  J  RT  In  10 


+  log 


RT  P° 


10 


N . 


(13) 


By  subtracting  Equation  12  from  13,  an  expression  for  relative  Henry’s  law 
constant  is  obtained: 


Mog10H 


ag:s;'  -ag:;^ 

j/w _ i/W 

2.303 RT 


(14) 


The  second  term  in  both  Equation  12  and  13  cancels  out  since  the  density 
of  pure  water  is  a  constant  at  any  specific  temperature.  Using  the  thermo¬ 
dynamic  path,  Equation  14  can  be  written  as: 


A Gt,  /z0  j)~  AG.,, ,  /z0  j) 

A  loo  H  —  Tr  W  J  _ Tr(vac)v 

~~  2.303 RT 


(15) 


The  absolute  Henry’s  law  constant  of  solute  z  is  then  calculated  from/s 
Henry’s  law  constant  using  equation: 

log10  H(  j)  =  A  log1Q  H  -  log1Q  H(  i) 


The  FEP  technique,  as  implemented  in  NAMD  simulation  engine  (Phillips 
et  al.  2005),  was  used  in  the  NPT  ensemble  for  computing  the  partition 
coefficients.  NAMD  uses  a  dual  topology  scheme  (Gao  et  al.  1989;  Pearl- 
man  1994),  where  both  the  initial  and  final  states  are  defined  concurrently. 
For  each  solute  of  interest,  three  FEP  simulations  were  performed  at 
298  K  and  1.013  bar;  one  for  the  water  phase,  one  for  the  water-saturated 
l-octanol  solution,  and  the  last  for  the  vacuum  phase. 

Simulations  were  also  run  at  308  and  318  K  to  investigate  the  temperature 
dependence  of  the  partition  coefficients.  The  mole  fraction  of  water  in  the 
octanol  phase  was  set  to  the  experimental  value  of  0.255  (Debolt  et  al. 
t995)-  FEP  was  carried  out  over  20  windows  where  the  starting  six  and  the 
ending  six  windows  were  unequally  spaced  with  very  small  increments  to 
improve  convergence  at  the  end  points.  This  methodology  is  known  to 
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avoid  the  end-point  catastrophe  (Beutler  et  al.  1994;  Pitera  et  al.  2002)  re¬ 
sulting  from  the  appearing  and  vanishing  atoms.  The  windows  between  0.1 
to  0.9  were  equally  spaced  at  0.1  increments.  A  non-bonded  cutoff  of  14  A 
and  a  timestep  of  1.0  fs  was  used.  The  Langevin  piston  Nose-Hoover 
method  (Martyna  et  al.  1994;  Feller  et  al.  1995)  was  used  to  control  pres¬ 
sure  and  temperature. 

The  Particle  Mesh  Ewald  (PME)  technique  (Essman  et  al.  1995)  was  used  to 
calculate  coulombic  interactions  in  all  MD  simulations.  Simulations  in  all 
phases  were  equilibrated  for  1  ns  before  free  energy  calculations  were  initi¬ 
ated.  For  the  calculations  in  vacuum,  an  isolated  hybrid  molecule  was  simu¬ 
lated  without  boundary  conditions  and  a  damping  coefficient  of  10  ps-1  for 
Langevin  temperature  control.  The  vacuum  run  was  carried  out  for  a  total  of 
2.4  ns  with  400  ps  of  equilibration  and  2  ns  of  sampling.  The  in  vacuo  simu¬ 
lations  require  relatively  longer  sampling  times  than  do  solvent  simulations. 
In  water  and  octanol  phases,  FEP  calculations  were  run  for  a  total  of  6  ns 
with  100  ps  of  equilibration  and  100  ps  of  sampling  for  each  window. 

Three  independent  simulation  trajectories  were  performed  in  each  phase 
and  the  values  averaged  for  the  net  free  energy  of  transfer,  which  is  used  to 
calculate  the  partition  coefficients.  Each  complete  FEP  simulation  re¬ 
quired  576  CPU  hours,  running  on  2.66  GHz  Intel  “Clovertown”  CPUs,  for 
simulations  of  solutes  in  water,  while  similar  calculations  performed  in 
water-saturated  octanol  required  960  CPU  hours  on  similar  hardware. 

4.2  Vapor-liquid  equilibria  and  vapor  pressure 

Gibbs-Duhem  integration  (Kofke  1993)  was  used  to  determine  the  phase 
coexistence  curve  (temperature  vs  density)  and  the  vapor  pressure.  With 
the  knowledge  of  an  initial  coexistence  point,  the  Clapeyron  equation  is 
integrated  to  provide  an  estimate  of  coexistence  points  at  other  tempera¬ 
tures.  The  Clapeyron  equation  is  given  by: 

dlnP  Ah 

-  = -  (17) 

dp  ]a  pPAv 

where: 

P  =  pressure 
P  =  l/kT 

Ah  =  difference  in  molar  enthalpies  of  the  coexisting  phases 

Av  =  difference  in  molar  volumes 
cr  indicates  that  the  derivative  is  taken  along  the  saturation  line. 
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The  method  allows  for  the  prediction  of  the  saturation  pressure  at  a  tem¬ 
perature  AT  away  from  the  known  coexistence  point  as  well.  Given  an  es¬ 
timate  of  the  saturation  pressure,  NPT  MD  simulations  are  performed 
simultaneously  for  both  liquid  and  vapor  phases  to  determine  the  coexist¬ 
ence  densities  and  heat  of  vaporization.  The  initial  coexistence  point  was 
determined  by  two  different  methods:  Grand  Canonical  Monte  Carlo 
(GCMC)  with  histogram  reweighting  technique  (Ferrenberg  et  al.  1988; 
Ferrenberg  et  al.  1989;  Potoff  et  al.  1998)  and  Performance  Verification 
Test  (PVT)  calculation  through  NPT  MD  simulation  near  the  critical  point. 
In  GCMC,  the  insertion  of  molecules  was  enhanced  through  multiple  first 
bead  insertions  and  the  application  of  the  coupled-decoupled 
configurational-bias  Monte  Carlo  method  (Martin  et  al.  1999).  The  ratios 
of  attempted  moves  were  set  to  60%  particle  insertions/deletions,  10% 
configurational-bias  regrowths,  15%  translations  and  15%  rotations. 

For  PVT  calculations,  isotherms  were  generated  at  different  pressures  near 
the  critical  point  and  densities  were  estimated.  One  isotherm,  where  liquid 
and  gas  coexist  at  a  specific  pressure,  is  chosen  as  the  initial  coexistence 
(P,  T)  condition.  For  Gibbs-Duhem  integration,  subsequent  gas  and  liquid 
simulations  starting  from  the  initial  coexistence  point  were  run  at  low 
temperatures  for  1  ns  each  with  300  ps  of  equilibration  and  700  ps  of 
sampling.  The  first  coexistence  simulation  was  carried  out  by  integrating 
the  Clapeyron  equation  with  trapezoidal  rule,  followed  by  two  simulations 
with  mid-point  predictor-corrector  method.  All  subsequent  simulations 
used  the  higher  order  Adams  predictor-corrector  integration  scheme.  A 
non-bonded  cutoff  of  14  A  without  tail  corrections  was  used  for  all  coexist¬ 
ence  simulations. 

4.3  Solid  phase  calculations 

4.3.1  Crystal  density  and  lattice  parameters 

Force  fields  were  validated  by  generating  lattice  parameters  and  crystal 
densities  and  comparing  them  to  the  experiment.  These  calculations  re¬ 
quire  knowledge  of  the  experimental  crystal  structures.  For  the  crystal 
density  and  lattice  parameter  calculations,  initial  crystal  structures  were 
taken  from  the  Cambridge  crystallographic  database  (Bruno  et  al.  2002) 
and  replicated  in  x,  y,  and  z  directions  to  create  a  supercell.  NPT  MD  simu¬ 
lations  were  run  at  zero  pressure  and  298  K.  The  system  was  initially  heat¬ 
ed  from  5  K  to  the  target  temperature  of  298  K  using  a  simulated  anneal¬ 
ing  technique.  The  temperature  and  pressure  control  methodology  is  the 
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same  as  MD  simulations  discussed  in  prior  sections.  The  system  was  equil¬ 
ibrated  for  l  ns,  where  first  250  ps  was  used  for  equilibration,  followed  by 
750  ps  of  time  averaging  for  the  cell  volume.  The  average  volume  was  then 
used  to  calculate  the  crystal  density. 

4.3.2  Melting  point 

A  solid-liquid  interface  method  based  on  the  work  of  Watt  et  al.  (2004) 
and  Morris  et  al.  (2002)  was  used  to  determine  the  melting  point.  A  solid- 
liquid  interfacial  system  was  prepared  as  follows: 

1.  In  the  original  supercell,  33%  of  the  molecules  were  constrained  to  fixed 
coordinates  and  the  rest  of  the  molecules  were  allowed  to  move. 

2.  A  few  molecules  were  permanently  removed  from  the  movable  region  to 
create  a  solid-liquid  interface. 

3.  The  structure  is  then  subjected  to  MD  simulations  in  the  NVT  ensemble 
around  1000  K  for  200  ps  to  create  liquid  regions  adjacent  to  the  fixed 
zone. 

4.  The  final  configuration  of  this  run  is  then  used  for  NPT  simulation  at  tem¬ 
peratures  close  to  the  experimental  melting  point. 

5.  Subsequent  to  this,  MD  simulations  in  the  NVE  ensemble  are  used  for 
equilibration  and  sampling  of  temperature  and  pressure. 

This  final  step  is  repeated  several  times  by  changing  the  volume  of  the  cell. 
The  resulting  temperatures  and  pressures  are  plotted  and  a  linear  regres¬ 
sion  fit  is  made.  The  temperature  corresponding  to  the  atmospheric  pres¬ 
sure  is  the  melting  point. 
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5  Results  and  Discussion 

5.1  Partition  coefficients 

Table  3  lists  the  net  free  energies  associated  with  each  transformation  in 
water,  water-saturated  octanol,  and  vacuum.  The  superscripts  denote  iter¬ 
ation  number.  The  values  from  each  iteration  are  averaged  to  calculate  the 
partition  coefficients. 

Figure  15  shows  the  computed  free  energies  with  respect  to  the  scaling  pa¬ 
rameter  1  in  each  phase  for  the  nitrobenzene  to  MNA  transformation,  re¬ 
spectively.  The  free  energies  in  the  plot  are  averages  from  three  iterations. 


Table  3.  Free  energies  predicted  in  water,  water-saturated  octanol  and  vacuum.  All  AG  are 

reported  in  kcal/mol. 


Transformation  (w) 

AG  liffw) 

AG2lr(w) 

AG3lr(w) 

Average 

Nitrobenzene  -  DNAN  (1) 

-16.68 

-16.45 

-16.66 

-16.60  +  0.12 

Nitrobenzene  -  MNA  (2) 

-1.60 

-1.70 

-1.82 

-1.71  ±  0.11 

Pyrazole  -  DNP  (3) 

83.75 

83.64 

84.24 

83.87  ±  0.31 

Pyrazole  -  NTO  (4) 

7.16 

6.56 

6.42 

6.71  ±  0.39 

Imidazole  -  MTNI  (5) 

-84.14 

-83.38 

-83.97 

-83.83  ±  0.39 

Nitrobenzene  -  TATB  (6) 

64.26 

65.54 

- 

64.90+0.90 

Transformation  (oct) 

AGlTr(oct) 

AG2Tr(oct) 

AG3Tr(oct) 

Average 

1 

-16.15 

-16.75 

-16.20 

-16.37  ±  0.33 

2 

-1.99 

-1.91 

-1.86 

-1.92  ±  0.06 

3 

85.16 

85.62 

84.91 

85.23  ±  0.36 

4 

9.07 

9.56 

9.62 

9.41  ±  0.30 

5 

-82.80 

-83.10 

-84.24 

-83.38  ±  0.75 

6 

61.19 

61.63 

- 

61.41+0.31 

Transformation  (vac) 

AGlTr(vac) 

AG2Tr(vac) 

AG3Tr(vac) 

Average 

1 

-11.46 

-11.43 

-11.41 

-11.43  ±  0.02 

2 

-0.52 

-0.52 

-0.52 

-0.52  ±  0 

3 

82.85 

82.86 

82.85 

82.85  +  0 

4 

13.50 

13.51 

13.51 

13.51  +  0 

5 

-81.97 

-82.11 

-81.74 

-81.94  +  0.18 

6 

51.88 

51.88 

- 

51.8810.0 
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Figure  15.  Free  energy  change  for  transformation  of  nitrobenzene  to  MNA  in 
water  (black),  octanol  (red),  and  vacuum  (green)  at  298K  and  1.013  bar. 


Table  4  lists  the  results  of  the  convergence  calculations  for  MNA,  where  i 
to  j  denote  the  forward  perturbation  and  j  to  i  the  reverse  perturbation. 
Each  entry  is  the  average  of  three  iterations  performed  in  each  phase.  The 
magnitude  of  the  incremental  free  energies  at  each  X  increment  and  the 
net  free  energy  change  for  the  forward  and  reverse  FEP  simulations  agree 
well  with  each  other  so  hysteresis  is  negligible  and  the  simulations  are 
considered  converged.  The  change  in  sign  is  due  to  the  difference  in  the 
direction  of  simulation.  The  free  energy  plots  and  convergence  calculations 
are  only  reported  for  the  nitrobenzene  to  MNA  transformations,  but  they 
are  representative  of  other  transformations. 


Table  4.  Computed  free  energies  (kcal/mol)  in  FEP  simulations  for  MNA. 


Water 

Octanol 

Vacuum 

U 

itoj 

jto  i 

itoj 

jto  i 

itoj 

jto  i 

0,  0.1 

2.46 

-2.11 

1.59 

-1.51 

-0.05 

0.05 

0.1,  0.2 

0.13 

-0.13 

-0.04 

0.05 

-0.05 

0.05 

0.2,  0.3 

0.15 

-0.15 

-0.25 

0.21 

-0.05 

0.05 

0.3,  0.4 

0.28 

-0.28 

-0.33 

0.32 

-0.05 

0.05 

0.4,  0.5 

0.39 

-0.38 

-0.36 

0.38 

-0.05 

0.05 

0.5,  0.6 

0.48 

-0.48 

-0.38 

0.44 

-0.05 

0.05 

0.6,  0.7 

0.58 

-0.57 

-0.44 

0.49 

-0.05 

0.05 

0.7,  0.8 

0.66 

-0.66 

-0.50 

0.52 

-0.05 

0.05 

0.8,  0.9 

-0.80 

0.77 

-0.57 

0.53 

-0.05 

0.05 

0.9,  1.0 

-0.94 

0.92 

-0.65 

0.68 

-0.05 

0.05 
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The  absolute  octanol-water  partition  coefficient  of  reference  solutes  nitro¬ 
benzene,  pyrazole,  and  imidazole,  and  Henry’s  law  constant  of  nitroben¬ 
zene  and  imidazole  were  taken  from  the  literature  (Sangster  1997;  Schultz 
et  al.  1982;  Hine  et  al.  1975;  SIDS  Report  2003).  Since  no  direct  Henry’s 
law  constant  has  been  reported  in  literature  for  pyrazole,  it  was  calculated 
from  experimental  vapor  pressure  and  solubility  of  pyrazole  at  298  K  by: 
H=  p/S  (18) 

where: 


p  =  vapor  pressure 
S  =  solubility. 

The  vapor  pressure  of  pyrazole  at  298  K  is  3.638  Pa  (Jimenez  et  al.  1987) 
and  solubility  in  water  at  298  K  is  19.4  mol/kg  of  water  (Wiley  1967).  Tables 
5  and  6,  respectively,  list  the  octanol-water  partition  coefficients  and  Hen¬ 
ry’s  law  constants  predicted  for  the  six  energetic  materials,  and  values  pre¬ 
dicted  using  COSMOtherm  by  Toghiani  et  al.  (2008),  EPI  Suite  (USEPA 
2009)  and  experimental  data  (Boddu  et  al.  2008;  Boddu  et  al.  2008). 


Table  5.  Octanol-water  partition  coefficient  (log  K0w). 


Molecule 

Simulation 

Exp 

EPI 

COSMO 

DNAN 

1.68 

1.61 

1.70 

1.92 

MNA 

2.00 

2.10 

2.01 

0.80 

DNP 

-0.97 

- 

-0.30 

0.37 

NTO 

-1.99 

- 

-1.56 

-1.19 

MTNI 

-0.40 

- 

0.05 

1.64 

TATB 

-1.86 

- 

-1.28 

4.74 

Table  6.  Henry's  Law  constants  (log  H). 


Molecule 

Simulation 

EPI 

Experiment 

DNAN 

-6.80 

-3.25 

-4.91 

MNA 

-3.88 

-3.60 

-6.17 

DNP 

-6.37 

-8.62 

- 

NTO 

-11.99 

-10.77 

- 

MTNI 

-9.24 

-9.69 

- 

TATB 

-12.56 

-14.45 

- 
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The  partition  coefficients  from  simulation  are  calculated  by  averaging  the 
forward  perturbation  results.  The  octanol-water  partition  coefficients  pre¬ 
dicted  by  FEP  simulations  are  within  ±  o.i  log  units  of  experiment  for  both 
DNAN  and  MNA.  While  EPI  Suite  values  also  predict  octanol-water  parti¬ 
tion  coefficients  in  good  agreement  with  the  experiment,  predictions  from 
COSMOtherm  have  unsigned  errors  of  0.12  and  1.3  log  units  for  DNAN 
and  MNA,  respectively.  Although  KOWWIN  (octanol-water  partition  coef¬ 
ficient  prediction  module  in  EPI)  was  developed  with  a  training  set  of 
about  2500  molecules  and  has  been  tested  on  a  dataset  of  10200  com¬ 
pounds,  it  might  give  poor  predictions  for  energetic  materials  since  the 
training  set  does  not  include  many  explosive  components.  COSMOtherm 
predicts  values  that  deviate  significantly  from  predictions  of  both  molecu¬ 
lar  simulations  and  the  EPI  Suite. 

The  Henry’s  law  constant  predicted  for  MNA  agrees  closely  with  the  exper¬ 
iment  while  for  DNAN,  it  is  under  predicted  significantly.  For  DNAN,  the 
source  of  error  is  unclear,  since  the  same  model  was  used  to  successfully 
predict  log  KoW  and  the  boiling  point  to  within  10%  of  experiment.  Alt¬ 
hough  the  relative  partitioning  between  octanol  and  water  was  predicted 
correctly,  it  is  possible  that  the  model  overpredicts  the  solubility  of  DNAN 
in  water,  leading  to  a  reduced  value  of  the  Henry’s  law  constant.  Investiga¬ 
tion  of  this  problem  is  ongoing. 

The  EPI  Suite  underpredicts  Henry’s  law  constants  of  both  DNAN  and 
MNA.  This  is  anticipated  because  HENRYWIN  (Henry’s  law  constant  pre¬ 
diction  module  of  EPI)  relies  on  a  much  smaller  calibration  set  of  just  345 
compounds  (Meylan  et  al.  1995),  therefore  the  predictive  capabilities  of 
the  EPI  Suite  in  this  respect  are  more  limited. 

Values  of  Henry’s  law  constants  indicate  the  volatility  of  the  compound. 
Compounds  with  Henry’s  law  constants  greater  than  io_5  atm.ms/mol  (log 
H  >  -3.39)  are  considered  highly  volatile  (Montgomery  2000).  None  of  the 
energetic  materials  fall  into  this  category  and  hence  partition  preferably 
into  the  aqueous  phase.  These  findings  are  further  illustrated  by  plotting 
the  predicted  partition  coefficients  in  the  multimedia-mobility  plot  pro¬ 
posed  by  Gillett  (1983).  The  predicted  partition  coefficients  are  located  in 
the  heavy  concern  area  D,  which  is  characterized  by  direct  effects  in  the 
water  column:  leaching  to  and  flow  through  groundwater  and  plant  root 
uptake.  The  compounds  are  not  predicted  to  bioaccumulate  or  induce  any 
atmospheric  problems. 
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5.2  Temperature  dependence  of  partition  coefficients 

The  temperature  dependence  of  the  octanol -water  partition  coefficients 
and  Henry’s  law  constants  were  analyzed  for  DNAN  and  MNA  by  obtain¬ 
ing  their  values  at  two  additional  temperatures  (308,  318)  apart  from 
298  K.  FEP  simulations  for  DNAN  and  MNA  were  run  at  308  and  318  Kin 
all  three  phases. 

Table  7  lists  octanol-water  partition  coefficients  and  Henry’s  law  constants 
predicted  at  different  temperatures  for  DNAN  and  MNA,  along  with  exper¬ 
imental  data  (Boddu  et  al.  2008;  Boddu  et  al.  2008).  Figure  16  shows  a 
plot  of  log  Kow  vs  l/T. 

The  octanol-water  partition  coefficients  decrease  with  increase  in  tem¬ 
perature  while  Henry’s  law  constants  increase  with  increase  in  tempera¬ 
ture.  A  reverse  trend  was  observed  for  the  MNA  experimentally,  where  the 
Henry’s  constant  decreased  with  increasing  temperature,  although  the  de¬ 
crease  is  small  and  the  statistical  error  in  the  data  is  unknown. 


Table  7.  Temperature  dependence  of  Partition  Coefficients  for  DNAN  and  MNA. 


DNAN 

MNA 

Temp 

log  Kow 

log  H 

log  Kow 

log  H 

(K) 

Sim 

Exp 

Sim 

Exp 

Sim 

Exp 

Sim 

Exp 

298 

1.68 

1.61 

-6.80 

-3.25 

2.00 

2.10 

-3.88 

-3.60 

308 

1.63 

1.54 

-6.56 

-3.24 

1.95 

1.98 

-3.83 

-3.64 

318 

1.54 

1.47 

-6.47 

-3.23 

1.92 

1.93 

-3.80 

-3.68 

Figure  16.  Octanol-water  partition  coefficient  as  a  function  of 
reciprocal  temperature  for  DNAN  (circle)  and  MNA  (square).  Filled 
symbols  correspond  to  experimental  values.  Solid  line  corresponds 
to  the  linear  regression  fit  to  simulation  data. 
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In  general,  experiments  have  shown  the  same  trend  as  the  simulation  data, 
i.e.,  Henry’s  law  constants  increase  with  temperature  increase  (Staudinger 
et  al.  1996).  Therefore  it  would  be  advisable  to  perform  additional  experi¬ 
ments  to  identify  the  source  of  the  unique  behavior  for  the  MNA  Henry’s 
constant  with  respect  to  temperature. 


The  data  in  the  plots  were  fit  to  van’t  Hoff  equation  (isochore),  which  gov¬ 
erns  the  variation  of  the  equilibrium  constant  with  temperature.  As  an 
equilibrium  constant,  log  Kow  can  be  expressed  as: 


l°9Ko 


-A H 

2.303 RT 


AS 

2.303 R 


(19) 


where: 

AH  =  enthalpy  of  water-octanol  partitioning 
AS  =  entropy  of  water-octanol  partitioning. 


Enthalpy  and  entropy  are  constant  over  the  temperature  range  studied. 

AH  and  AS  are  determined  from  the  linear  regression  fit  to  the  log  KoW  vs 
l/T  plot.  The  Gibbs  free  energy  of  partitioning  (AG)  at  a  specific  tempera¬ 
ture  is  determined  from  Equation  5.  Table  8  lists  the  enthalpy,  entropy, 
and  Gibbs  free  energy  of  partitioning  for  DNAN  and  MNA  between  octanol 
and  water,  along  with  the  experimental  values.  For  both  DNAN  and  MNA, 
transfer  from  water  to  octanol  is  exothermic  and  enthalpy  driven,  which  is 
evident  from  the  negative  values  of  AH.  The  temperature  dependence  of 
the  Henry’s  law  constant  is  described  by  the  equation: 


log  H 


-AH 

_ V 

2.303 RT 


-AS 

_ V__ 

2.303 R 


(20) 


where: 

AHV  =  Enthalpy  of  volatilization 
ASV  =  Entropy  of  volatilization. 

Figure  17  shows  (and  Table  9  lists)  the  enthalpy  and  entropy  of  volatiliza¬ 
tion  obtained  from  the  linear  regression  fit  to  log  H  vs  l/T  plot. 
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Figure  17.  Henry’s  law  constant  as  a  function  of  reciprocal  temperature  for 
DNAN  (circle)  and  MNA  (square).  Filled  symbols  correspond  to  experimental 
values.  Solid  line  corresponds  to  the  linear  regression  fit  to  simulation  data. 


Table  8.  Free  energy,  enthalpy,  and  entropy  of  water-octanol  partitioning. 


Property 

DNAN 

MNA 

Sim 

Exp 

Sim 

Exp 

AG298  K(kJ/mol) 

-9.58 

-9.22 

-11.41 

-11.95 

AH  (kJ/mol) 

-12.65 

-12.70 

-7.27 

-15.06 

AS  (J/mol/K) 

-10.27 

-11.68 

13.83 

-10.44 

Table  9.  Enthalpy  and  entropy  of  water-air  partitioning. 


Property 

DNAN 

MNA 

Sim 

Exp 

Sim 

Exp 

AHv  (kJ/mol) 

30.06 

2.15 

3.15 

-6.62 

ASv  (J/mol/K) 

-28.86 

-55.18 

-21.63 

-91.30 

The  positive  enthalpy  change  indicates  that  transfer  from  water  to  gaseous 
state  is  an  endothermic  process.  Negative  entropy  of  transfer  and  a  posi¬ 
tive  enthalpy  term  suggest  that  volatilization  is  neither  enthalpy  nor  en¬ 
tropy  driven  (the  process  is  not  spontaneous)  and  the  compounds  have 
strong  interactions  in  aqueous  solution. 

Vapor-Liquid  Equilibria  and  Vapor  Pressure 

The  force  field  developed  for  these  compounds  can  also  be  used  to  com¬ 
pute  other  properties  like  critical  parameters,  boiling  points,  vapor  pres¬ 
sure,  heats  of  vaporization,  and  acentric  factor.  The  vapor-liquid  coexist¬ 
ence  curves  and  critical  parameters  are  useful  in  equation  of  state 
modeling  of  these  compounds.  Moreover,  these  properties  can  be  used  in 
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the  development  of  QSPR/QSAR  models,  in  contrast  to  boiling  points  and 
critical  parameters  derived  from  empirical  correlations,  to  improve  their 
predictive  capability.  Figures  18  and  19,  respectively,  show  vapor-liquid 
coexistence  curves  and  vapor  pressure  plots  for  DNAN  and  MNA.  The 
phase  diagrams  for  DNAN  and  MNA  should  be  considered  hypothetical, 
since  these  compounds  are  known  to  decompose  at  temperatures  near 
their  normal  boiling  points. 


Figure  18.  Vapor-liquid  coexistence  curves  for  DNAN  (circle)  and  MNA  (square).  Line  is  a  fit  of 
simulation  data  to  scaling  laws.  Filled  symbols  correspond  to  predicted  critical  points. 


Figure  19.  Clausius-Clapeyron  plot  for  DNAN  (circle)  and  MNA  (square). 


1/T  (K1) 
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Critical  temperatures  and  densities  were  computed  by  fitting  the  saturated 
liquid  and  vapor  densities  to  the  density  scaling  law  for  critical  tempera¬ 
ture  (Rowlinson  et  al.  1982): 

p«t-p  mr=B<T-T,r  (2D 

and  the  law  of  rectilinear  diameters  (Rowlinson  et  al.  1982): 

Pr  +P 

bq  vap  =pc  +  A(T-T) 

2  ‘  ‘  (22) 

where: 

/?  =  critical  exponent  for  Ising-type  fluids  in  3  dimensions  (0.325) 
A,B  =  constants  fit  to  simulation  data. 

Table  10  lists  the  critical  parameters,  boiling  point,  and  acentric  factor, 
along  with  values  predicted  through  group  contribution  method  (Stein  et 
al.  1994)  by  Toghiani  et  al.  (Toghiani  et  al.  2008).  The  experimental  boil¬ 
ing  point  of  DNAN  at  12  mm  Hg  is  479  K  (CRC  Handbook  of  Chemistry  & 
Physics,  2008-2009).  The  vapor  pressure  data  from  simulation  were  ex¬ 
trapolated  using  the  Clausius-Clayperon  equation  to  12  mm  Hg  (0.016 
bar).  The  temperature  corresponding  to  this  pressure  is  461.04  K,  about 
3.7%  lower  than  the  experiment.  The  difference  between  the  values  pre¬ 
dicted  by  simulation  and  group  contribution  is  more  pronounced  for 
DNAN  than  MNA,  which  may  be  due  to  the  proximity  of  the  ortho-nitro 
group  and  the  methoxy  group  in  DNAN.  In  the  force  fields  developed  in 
this  work,  the  proximity  of  other  functional  groups  and  synergistic  effects 
are  taken  into  account  through  the  partial  charge  distributions,  which  are 
derived  from  electrostatic  potential  energy  surfaces  determined  from  ab 
initio  calculations.  For  MNA,  the  nitro  and  the  amine  groups  are  far 
enough  apart  that  synergy  effects  are  expected  to  be  negligible. 

Table  10.  Critical  parameters,  boiling  point,  and  acentric  factor  for  DNAN  and  MNA. 


Molecule 

To  (K) 

pc  (kg/m3) 

Pc  (bar) 

Tb  (K) 

00 

DNANa 

885.42 

410.20 

37.36 

620.82 

1.54 

DNANb 

806 

- 

39.90 

588 

0.85 

MNAa 

770.75 

324.50 

37.70 

522.76 

1.41 

MNAb 

748 

- 

41.70 

527 

0.65 

a  This  work 

b  Group  Contribution 
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Heats  of  vaporization  were  calculated  for  each  molecule  as  a  function  of 
temperature  using  Gibbs-Duhem  integration  data  collected  for  the  vapor- 
liquid  equilibria  calculations  and  Equation  23: 

^-U-U^P^-V,)  (23) 

where: 

U  =  internal  energy  per  mol 
V  =  molar  volume. 

The  subscripts  v  and  l  refer  to  the  vapor  and  liquid  phases,  respectively. 
Figure  20  shows  the  results  of  these  calculations. 

5.3  Solid  phase  calculations 

Tables  11  and  12  list  the  crystal  lattice  parameters  and  density  determined 
for  all  the  six  compounds. 

The  lattice  parameters  and  crystal  density  predicted  for  DNAN,  MNA  and 
NTO  are  in  good  agreement  with  the  experiment.  For  MTNI,  simulation 
underpredicts  the  c  dimension  and  slightly  over  predicts  the  crystal  densi¬ 
ty.  Since  experimental  crystal  structure  is  not  available  for  DNP,  solid 
phase  calculations  were  not  performed  for  it. 

Melting  Point 

The  melting  point  of  NTO  was  determined  using  solid-liquid  interface 
method  as  discussed  earlier.  Figure  21  shows  the  initial  configuration  used 
for  NVE  ensemble.  Figure  22  shows  the  temperature-pressure  plot. 

Although  pressure  and  temperature  do  not  have  a  linear  relation,  the  small 
range  of  temperatures  (530-560  K)  covered  allows  us  to  assume  a  linear 
dependence.  The  melting  point  predicted  is  538.69  K,  which  is  in  excellent 
agreement  with  the  experimental  value  of  539.35  K  (Liu  et  al.  1995). 
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Figure  20.  Heat  of  vaporization  for  DNAN  (circle)  an  MNA  (square) 
predicted  from  NPT  MD  simulations. 


900 


Table  11.  Crystal  parameters  and  density  for  DNAN,  MNA,  and  NTO. 


Parameters 

DNAN  (Monoclinic) 

MNA  (Monoclinic) 

NTO  (Triclinic) 

Sim 

Exp1 

Sim 

Exp2 

Sim 

Exp3 

a  (A) 

9.15 

8.77 

9.78 

10.07 

5.21 

5.12 

b  (A) 

12.23 

12.64 

7.02 

6.93 

10.50 

10.30 

c  (A) 

15.63 

15.42 

11.07 

10.81 

18.32 

17.9 

a 

90 

90 

90 

90 

106.58 

106.7 

p 

81.64 

81.89 

101.32 

101.95 

97.79 

97.7 

Y 

90 

90 

90 

90 

90.11 

90.2 

p  (g/cm3) 

1.52 

1.56 

1.36 

1.36 

1.81 

1.92 

1  (Nyburg  et  al.  1987) 

2  (Schaefer  et  al.  1988) 

3  (Bolotina  et  al.  2005) 

Table  12.  Crystal  parameters  and  density  of  MTNI  and  TATB. 


MTNI  (Orthorhombic) 

TATB  (Triclinic) 

Parameters 

Sim 

Exp1 

Sim 

Exp2 

a  (A) 

8.51 

8.61 

8.89 

9.01 

b  (A) 

17.70 

17.71 

8.91 

9.02 

c  (A) 

9.89 

10.68 

6.64 

6.81 

a 

90 

90 

108.77 

108.59 

p 

90 

90 

91.82 

91.82 

Y 

90 

90 

119.95 

119.97 

P  (g/cm3) 

1.92 

1.76 

2.03 

1.93 

1  (Cho  et  al.  2002) 

2  (Cady  et  al.  1965) 
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Figure  21.  Snapshot  of  initial  configuration  used  for  NVE  simulations  of 
pressure  and  a  linear  regression  fit  is  made.  The  temperature  corresponding 
to  the  atmospheric  pressure  is  the  melting  point. 


Figure  22.  Temperature-pressure  plot  for  the  coexistence  system. 


ERDC/CERL  TR-10-26 


33 


6  Conclusion 

This  work  has  demonstrated  the  potential  of  atomistic  computer  simula¬ 
tions  for  the  prediction  of  partitioning  and  physical  property  prediction  for 
energetic  materials.  Force  fields  were  developed  for  six  energetic  materials 
(DNAN,  MNA,  DNP,  NTO,  MTNI  and  TATB)  and  the  predicted  thermo¬ 
physical  properties  were  found  to  be  in  close  agreement  (5-10%  in  most 
cases)  with  the  scarce  experimental  data  available.  Based  on  the  predicted 
octanol -water  and  Henry’s  law  constants,  with  the  exception  of  TATB,  all 
compounds  studied  in  this  work  are  predicted  to  be  problematic  with  re¬ 
spect  to  groundwater  contamination. 

In  addition  to  the  properties  calculated  in  this  report,  the  generalized, 
transferable  force  fields  for  energetic  materials  presented  here  may  be 
used  to  investigate  the  interactions  of  energetic  materials  in  a  wide  variety 
of  complex  systems,  including  their  diffusion  and  transport  in  the  envi¬ 
ronment. 
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Acronyms  and  Abbreviations 


Term 

Definition 

ARDEC 

US  Armament  Research,  Development,  and  Engineering  Center 

CERL 

Construction  Engineering  Research  Laboratory 

CHELPG 

CHarges  from  Electrostatic  Potentials  using  a  Grid  based  method 

COSMO 

Conductor  Screening  Module 

CPU 

central  processing  unit 

DC 

District  of  Columbia 

DNAN 

2.4-dinitroanisole 

DNP 

3,4-Dinitropyrazole 

EH 

Explicit  Hydrogen 

EPI 

Estimation  Program  Interface 

ERDC 

Engineer  Research  and  Development  Center 

FEP 

Free  Energy  Perturbation 

GB/SA 

Generalized  Born/Surface  Area 

GCMC 

Grand  Canonical  Monte  Carlo 

HF 

Hartree-Fock 

IM 

Insensitive  Munitions 

MC 

Monte  Carlo 

MD 

Molecular  Dynamics 

MNA 

n-Methyl-p-nitroaniline 

MTNI 

l-Methyl-2,4,5-trinitroimidazole 

NTO 

3-Nitro-l,2,4-triazol-5-one 

PI 

Principal  Investigator 

PME 

Particle  Mesh  Ewald 

PVT 

Performance  Verification  Test 

QCISD 

Quadratic  Configuration  Interaction  Singles  Doubles 

QSAR 

Quantitative  Structure-Activity  Relationship 

QSPR 

Quantitative  Structure-Property  Relationship 

SF 

Standard  Form 

SIDS 

Screening  Information  Dataset 

SM 

Solvation  Model 

TATB 

l,3,5-Triamino-2,4,6-trinitrobenzene 

TR 

Technical  Report 

UA 

United -Atom 

URL 

Universal  Resource  Locator 

US 

United  States 

USE  PA 

US  Environmental  Protection  Agency 

WWW 

World  Wide  Web 
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Appendix  A:  Lennard-Jones  Parameters  and 
Partial  Charges 


Table  Al.  Lennard-Jones  parameters  for  DNAN 
(UA). 


Site 

Molecule 

a  (A) 

e(K) 

q(e) 

C<x 

Ether 

4.50 

15 

0.150 

Ca 

Nitro 

4.50 

15 

0.112 

CH 

Benzene 

3.74 

48 

0.00 

CHs 

Ether 

3.75 

98 

0.252 

0 

Ether 

2.80 

55 

-0.402 

N 

Nitro 

3.31 

40 

0.768 

0 

Nitro 

2.90 

80 

-0.440 

Table  A2.  Lennard-Jones  parameters  for  MNA 
(UA). 


Site 

Molecule 

a  (A) 

e(K) 

q(e) 

Nitro 

4.50 

15 

0.131 

Amine 

4.50 

15 

0.187 

CH 

Benzene 

3.74 

48 

0.00 

CHs 

Amine 

3.75 

98 

0.234 

N 

Nitro 

3.31 

40 

0.711 

N 

Amine 

3.52 

58 

-0.730 

0 

Nitro 

2.9 

80 

-0.449 

H 

Amine 

0 

0 

0.365 

Table  A3.  Lennard-Jones  parameters  for  DNAN 
(EH). 


Site 

Molecule 

a  (A) 

e(K) 

q(e) 

Ca 

Ether 

3.60 

30.7 

0.150 

c„ 

Nitro 

3.60 

30.7 

0.090,  0.142 

C-(H) 

Ring 

3.60 

30.7 

-0.165,  -0.165,  -0.189 

H-C 

Ring 

2.36 

25.45 

0.165,  0.165,  0.189 

C 

Methyl 

3.6 

47 

0.132 

H 

Methyl 

2.5 

10 

0.041 

0 

Ether 

2.80 

55 

-0.407 

N 

Nitro 

2.90 

30 

0.774,  0.723 

0 

Nitro 

2.70 

42 

-0.432' 
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Table  A4.  Lennard-Jones  parameters  for  MNA  (EH). 
Site  Molecule  a  (A)  e  (K)  q(e) 


Site 

Molecule 

a  (A) 

e(K) 

q(e) 

N 

Nitro 

2.90 

30 

0.702 

0 

Nitro 

2.70 

42 

-0.414 

N 

Sp2 

3.20 

57 

-0.396 

C 

Nitro 

3.60 

30.7 

0.354 

N 

Amide 

3.40 

141 

-0.023 

H 

Amide 

0.50 

12 

0.321 

C 

Sp2 

3.60 

30.7 

-0.309 

H 

Attached  to  ring  carbon 

2.36 

25.45 

0.206 

Table  A6.  Lennard-Jones  parameters  for  NTO. 


Site 

Molecule 

ct(A) 

6(K) 

q(e) 

0 

Nitro 

2.70 

42 

-0.416 

N 

Nitro 

2.90 

30 

0.722 

C 

Sp2 

3.60 

30.7 

0.408 

N 

Sp2 

3.20 

57 

-0.387 

N 

Amidel(attached  to  sp2  nitrogen) 

3.40 

141 

-0.187 

H 

Amidel 

0.50 

12 

0.315 

C 

Carbonyl 

3.60 

30.7 

0.689 

0 

Carbonyl 

3.05 

79 

-0.601 

N 

Amide2  (attached  to  sp2  carbon) 

3.40 

141 

-0.476 

H 

Amide2 

0.50 

12 

0.349 
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Table  A7.  Lennard-Jones  parameters  for  MTNI. 


Site 

Molecule 

<j(A) 

e(K) 

q(e) 

N 

Nitro 

2.90 

30 

0.742 

0 

Nitro 

2.70 

42 

-0.404 

N 

Sp2 

3.20 

57 

-0.529 

C 

Methyl 

3.75 

98 

0.236 

N 

Ring 

3.40 

141 

-0.047 

C 

Nitro 

3.60 

30.7 

0.403,  -0.199,  0.334 

Table  A8.  Lennard-Jones  parameters  for  TATB. 


Site 

Molecule 

a  (A) 

s(K) 

q(e) 

c„ 

Nitro 

3.60 

30.7 

0.061 

c„ 

Amine 

3.60 

30.7 

0.076 

N 

Nitro 

2.90 

30 

1.131 

N 

Amine 

3.26 

160 

-1.110 

0 

Nitro 

2.70 

42 

-0.548 

H 

Amine 

0.50 

12 

0.519 

Table  A9.  Bond  parameters  for  DNAN. 


Bond 

Molecule 

Equilibrium  Bond 
Length  (A) 

Force  constant 
(kcal/mol) 

0-N 

nitro 

1.22 

866.45 

N-C 

nitro,  aromatic 

1.49 

363.08 

C-C 

aromatic,  aromatic 

1.4 

528.27 

C-0 

aromatic,  ether 

1.41 

480.35 

o-c 

ether,  methyl 

1.41 

289.56 

Table  A10.  Bond  parameters  for  MNA. 


Bond 

Molecule 

Equilibrium  Bond  Length  (A) 

Force  constant  (kcal/mol) 

0-N 

nitro 

1.22 

872.54 

N-C 

nitro,  aromatic 

1.49 

361.61 

C-C 

aromatic,  aromatic 

1.4 

529.35 

C-N 

aromatic,  amine 

1.35 

528.94 

N-C 

amine,  methyl 

1.44 

413.41 

N-H 

amine 

0.99 

614.35 
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Table  All.  Bond  parameters  for  DNP. 


Bond 

Equilibrium 
Bond  Length  (A) 

Force  constant 
(kcal/mol) 

C-C 

1.40 

485.9 

C=N 

1.30 

646.2 

C-Nitro 

1.43 

409.6 

N-C 

1.34 

566.2 

N-H 

0.99 

601.8 

N-0 

1.18 

956.5 

C=C 

1.35 

636.4 

N-N 

1.30 

557.8 

Table  A12.  Bond  parameters  for  NTO. 


Bond 

Molecule 

Equilibrium  Bond 
Length  (A) 

Force  constant 
(kcal/mol) 

H-N 

Amide 

0.99 

611.64 

N-C 

Amide,  Carbonyl 

1.37 

428.83 

N-N 

Amide,  Sp2 

1.35 

435.90 

N=C 

Sp2 

1.25 

932.04 

C-N 

Sp2,  Amide 

1.35 

459.76 

C=0 

Carbonyl 

1.19 

1061.00 

C-N 

Sp2,  Nitro 

1.44 

377.16 

N-0 

Nitro 

1.18 

1041.04 

Table  A13.  Bond  parameters  for  MTNI. 


Bond 

Molecule 

Equilibrium  Bond 
Length  (A) 

Force  constant 
(kcal/mol) 

N-0 

nitro 

1.17 

1007 

N-C 

nitro 

1.44 

388.30 

C=C 

ring 

1.35 

621.20 

C=N 

ring 

1.27 

723.70 

C-N 

ring 

1.34 

535.40 

N-C 

sp2  nitrogen 

1.34 

535.4 

N-C 

sp2  carbon 

1.33 

529.4 
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Table  A15.  Bending  parameters  for  DNAN. 


Angle 

Molecule 

Equilibrium  Bond 
Angle  (degree) 

Force  constant 
(kcal/mol) 

C-O-C 

ether,  aromatic 

122 

97.94 

O-N-O 

Nitro 

125 

181.13 

O-N-C 

nitro,  aromatic 

117.5 

167.89 

C-C-C 

aromatic,  aromatic,  aromatic 

120 

189.7 

O-C-C 

ether,  aromatic,  aromatic 

125 

138.72 

N-C-C 

nitro,  aromatic,  aromatic 

120 

154.79 

Table  A16.  Bending  parameters  for  MNA. 


Angle 

Molecule 

Equilibrium  Bond 
Angle  (degree) 

Force  constant 
(kcal/mol) 

H-N-C 

amine 

117.77 

72.9 

O-N-O 

nitro 

125 

181.1 

O-N-C 

nitro,  aromatic 

117.5 

167.9 

C-C-C 

aromatic,  aromatic,  aromatic 

120 

189.4 

N-C-C 

amine,  aromatic,  aromatic 

120 

145.4 

N-C-C 

nitro,  aromatic,  aromatic 

120 

154.8 
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Table  A17.  Bending  parameters  for  DNP. 


Bond 

Equilibrium  Bond 
Angle  (degree) 

Force  constant 
(kcal/mol) 

H-N-C 

127.24 

75.22 

N-C=C 

108.77 

296 

N=C-Nitro 

120.82 

130.6 

O-N-O 

126.7 

182.2 

N-N-C 

101.5 

317 

c=c-c 

101.21 

317 

C-C=N 

113.43 

290.5 

C=C-Nitro 

130.73 

107.8 

C-C-Nitro 

125.73 

123.2 

C=N-N 

111.97 

322.8 

C-N-0 

117.16 

144.8 

H-N-N 

121.14 

80.65 

N-C-Nitro 

120.49 

122 

Table  A18.  Bending  parameters  for  NTO. 


Bond 

Molecule 

Equilibrium  Bond 
Angle  (degree) 

Force  constant 
(kcal/mol) 

C-N-C 

sp2,  amidel,  carbonyl 

106.88 

2514.45 

H-N-C 

amide2,  carbonyl 

125.93 

83.645 

H-N-N 

amide2,  sp2 

120.44 

87.99 

H-N-C 

amidel,  carbonyl 

125.94 

80.57 

H-N-C 

amidel,  sp2 

127.18 

73.36 

N-C-N 

Nitro,sp2,  amidel 

121.55 

149.87 

N-N-C 

amide2,  sp2, 

103.74 

2401.53 

N-C-N 

sp2,  sp2,  Nitro 

124.47 

149.89 

N-C-N 

amidel,  carbonyl, amide2 

101.77 

1347.10 

N-C-0 

amidel,  carbonyl 

129.18 

136.36 

0-C-N 

carbonyl,  amide2 

129.05 

231.29 

O-N-O 

Nitro 

127.18 

182.20 

0-N-C 

Nitro,  sp2 

118.10 

231.29 
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Table  A19.  Bending  parameters  for  MTNI. 


Bond 

Molecule 

Equilibrium  Bond 

Angle  (degree) 

Force  constant 
(kcal/mol) 

O-N-O 

Nitro 

126.70 

182.20 

O-N-C 

Nitro 

117.12,  116.47,  118.01 

148.50,  116.50,  140.20 

N-C=C 

Nitro, Csp2 

131.50 

110.80 

N-C-N 

Ring,  Nitro 

123.13 

160.40 

C=C-N 

Ring,  Amide 

107.37 

296.70 

C-N-C 

Ring 

103.50 

337.70 

C-N-C 

Ring,  Methyl 

130.6 

130.60 

N-C=N 

Ring,  sp2 

114.46 

283.70 

Table  A20.  Bending  parameters  for  TATB. 


Angle 

Molecule 

Equilibrium  Bond 
Angle  (degree) 

Force  constant 
(kcal/mol) 

H-N-C 

amine 

119.80 

73.86 

O-N-O 

nitro 

125.00 

181.13 

O-N-C 

nitro,  aromatic 

117.5 

167.91 

C-C-C 

aromatic,  aromatic,  aromatic 

120 

189.41 

N-C-C 

amine,  aromatic,  aromatic 

120 

145.43 

N-C-C 

nitro,  aromatic,  aromatic 

120 

154.81 

Table  A21.  Torsional  parameters  for  DNAN. 


Dihedral 

Molecule 

n 

Phase  angle 
(degree) 

0  [kcal/mol] 

C-C-C-C 

aromatic 

2 

180 

15.230 

0-N-C-C  (ortho) 

nitro,  aromatic 

1,  2,  3,  4 

0,  0,  0,  0 

0.065,  -0.202,  0.085,  0.571 

0-N-C-C  (para) 

nitro,  aromatic 

1,  2 

180,  180 

-0.136,  4.351 

C-O-C-C 

ether,  aromatic 

1,  2 

180,  180 

0.663,  1.467 

Table  A22.  Torsional  parameters  for  MNA. 


Dihedral 

Molecule 

n 

Phase  angle 
(degree) 

G  [kcal/mol] 

C-C-C-C 

Aromatic 

2 

180 

15.230 

0-N-C-C  (ortho) 

nitro,  aromatic 

1,  2 

180,  180 

-0.136,  4.351 

C-N-C-C 

amine, 

aromatic 

2,  4 

180,  180 

-0.308,  3.003 
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Table  A23.  Torsional  parameters  for  DNP. 


Dihedral 

Molecule 

n 

Phase  angle 
(degree) 

Ci  [kcal/mol] 

C=C-N-N,  C-N-N=C 

Ring 

1 

180 

111.600 

N-N=C-C,  N=C-C=C 

Ring 

1 

0 

134.400 

C-C=C-N 

Ring 

1 

180 

144.00 

O-N-C-C 

Nitro 

1,2 

180 

-0.082,  3.29 

Table  A24.  Torsional  parameters  for  NTO. 


Dihedral 

Molecule 

N 

Phase  angle 
(degree) 

Ci  [kcal/mol] 

C-C=N-N 

ring 

1 

180 

50.62 

C=N-N-C,  N-C-N-C,  C-N-C=N 

ring 

1 

180 

69.83 

N-N-C-N 

ring 

1 

180 

104.60 

O-N-C-C 

nitro 

1,2 

180,  180 

-0.082,  3.29 

Table  A25.  Torsional  parameters  for  MTNI. 


Dihedral 

Molecule 

n 

Phase  angle 
(degree) 

Ci  [kcal/mol] 

N=C-N-C 

Ring 

1 

180 

123.40 

C-N-C=C 

Ring 

1 

180 

143.80 

N-C=C-N 

Ring 

1 

180 

125.30 

C=C-N=C 

Ring 

1 

180 

134.80 

C-N=C-N 

Ring 

1 

180 

126.20 

0-N-C-N 

Nitro 

1,  2 

180,  180 

-0.059,  1.218 

0-N-C=C 

Nitro 

1,  2 

0,  0 

0.065,  0.584 

Table  A26.  Torsional  parameters  for  TATB. 


Dihedral 

Molecule 

n 

Phase  angle  (degree) 

Ci  [kcal/mol] 

C-C-C-C 

ring 

2 

180 

15.230 

O-N-C-C 

nitro 

1,  2,  3,  4 

180,  180,  180,  180 

0.023,  4.755,  0.017,-1.014 

H-N-C-C 

amine 

1,  2,  3,  4 

0,  180,  0,  180 

0.023,  3.015,  -0.451, 

0.119 
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gest  that  empirical  force  fields,  combined  with  molecular  dynamics  simulations,  provide  an  accurate  methodology  for  predicting  rele¬ 
vant  descriptors  of  environmental  fate  for  energetic  materials. 
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